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In this investigation we revisit the concept of "effective free surfaces" arising in 
the solution of the time-averaged fluid dynamics equations in the presence of free 
boundaries. This work is motivated by applications of the optimization and opti- 
mal control theory to problems involving free surfaces, where the time-dependent 
formulations lead to many technical difficulties which are however alleviated when 
steady governing equations are used instead. By introducing a number of pre- 
cisely stated assumptions, we develop and validate an approach in which the 
interface between the different phases, understood in the time-averaged sense, is 
sharp. In the proposed formulation the terms representing the fluctuations of the 
free boundaries and of the hydrodynamic quantities appear as boundary condi- 
tions on the effective surface and require suitable closure models. As a simple 
model problem we consider impingement of free-falling droplets onto a fluid in 
a pool with a free surface, and a simple algebraic closure model is proposed for 
this system. The resulting averaged equations are of the free-boundary type and 
an efficient computational approach based on shape optimization formulation is 
developed for their solution. The computed effective surfaces exhibit consistent 
dependence on the problem parameters and compare favorably with the results 
obtained when the data from the actual time-dependent problem is used in lieu 
of the closure model. 
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I. INTRODUCTION 



The goal of this work is to investigate the concept of "effective free surfaces" which 
are defined here as stationary interfaces corresponding to the time-averaged balance 
of mass, momentum and, if applicable, energy in a time-dependent flow with free sur- 
faces. In other words, given an unsteady two-phase flow with fluctuating boundaries, 
the effective free surface represents the boundary between the two phases in the cor- 
responding mean flow which satisfies the time-averaged form of the original system of 
governing equations subject to a number of modeling assumptions. The motivation for 
this work comes from the field of flow control 1 where many emerging applications involve 
control and optimization of free-surface phenomena. The particular applications under- 
lying this research concern optimization of complex thermo-fluid phenomena occurring 
in liquid metals during welding, see Volkov et a/. 2 . While the mathematical founda- 
tions for the optimal control of time-dependent free-boundary problems are relatively 
well understood 3 , such approaches tend to result in computational problems of signif- 
icant complexity even for simple models 4 . The main difficulty arising when methods 
of the optimal control, or more broadly, the calculus of variations are applied to such 
problems is that some of the optimality conditions have the form of partial differential 
equations (PDEs) defined on interfaces which evolve with time. Needless to say, such 
problems tend to be quite hard to solve for non-trivial applications. On the other hand, 
this framework becomes much more tractable when time-independent free-boundary 
problems are considered instead 5 . Moreover, on a more practical level, fluid flows with 
free surfaces may generate "subgrid-scale" features which are particularly difficult to 
compute, and it is therefore desirable to account for their effect in the average balance 
of mass and momentum in a systematic manner. In this paper we propose and test 
a simple mathematical model, in the form of a system of coupled PDEs of the free- 
boundary type, representing the time-averaged conservation of mass and momentum in 
a given time-dependent problem with free surfaces. While such averaging approaches 
are well-established in the study of turbulent flows in domains with fixed boundaries, 
giving rise to the well-known Reynolds-Averaged Navier-Stokes (RANS) equations, see, 
e.g., Pope 6 , the additional complication in the present problem is that one also needs to 
take into account the effect of the fluctuations of the location of the free surfaces on the 
average mass and momentum balance. Our approach to this problem relies on a number 
of simplifying assumptions which are all clearly identified. In the spirit of the "closure 
problem" arising in turbulence modeling, see Ref. 6, in order to close the resulting sys- 
tem of equations one needs to express average products of fluctuating quantities in terms 
of average quantities. However, in contrast to the classical closure problem where the 
Reynolds stresses are modeled with terms defined in the bulk of the fluid, in the present 
problem, subject to certain assumptions, such closure terms will appear in the bound- 
ary conditions defined on the effective free boundary. We will also discuss some very 
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simple strategies for constructing such closures. The question of ensemble-averaged, or 
time-averaged, description of flows with interfaces has received some attention in the 
literature and we mention here the work of Dopazo 17 , Hong & Walker 18 and Brocchini & 
Peregrine 19 ' 20 which also contains references to a number of earlier attempts. These prob- 
lems were recently revisited in the context of the derivation and validation of suitable 
models for multiphase Reynolds-averaged Navier-Stokes (RANS) equations 21 and Large 
Eddy Simulation (LES) 22-24 . We also mention the recent investigation by Waclawczyk 
& Oberlack 25 where a number of closure strategies were proposed for this type of flows. 
Finally, we add that the related question of homogenization of free-boundary problems 
is an emerging topic in the mathematical analysis of PDEs, see, e.g., Schweizer 26 . A 
detailed description of various computational methods applied to multiphase flows can 
be found in the monograph by Prosperetti & Tryggvason 27 . As compared to these earlier 
studies, novel aspects of the present investigation are that, first, we want to compute 
steady-state solutions, which is motivated by the optimization applications mentioned 
above, and secondly, we want our averaged flows to feature sharp effective surfaces, so 
that the free-boundary property of the original problem is preserved in its averaged 
version. In contrast, we note that the formulations developed in Refs. 20 and 25 lead to 
interfaces, referred to as "surface layers" , characterized by finite thickness. We also wish 
to highlight that although Brocchini & Peregrine 20 derived averaged equations taking 
into account the fluctuations of the free boundaries and also proposed a simple closure 
model, to the best of our knowledge, there have been no attempts to actually compute 
such solutions for nontrivial problems which is one of the contributions of the present 
work. 

The resulting system of PDEs represents the averaged balance of mass and momentum 
which has the form of a steady-state free-boundary problem. Since such problems tend 
to be difficult to solve numerically, we also propose a solution approach based on shape 
optimization which is well adapted to the numerical solution of this class of problems. In 
order to test our approach we choose a very simple model problem which, while allowing 
us to focus on certain methodological aspects, still captures some essential features of the 
motivating application, namely, the transfer of mass and momentum to the weld pool via 
droplets, see Figure 1. This model describes the two-dimensional (2D), time-periodic 
impingement of droplets on the free surface of the fluid in a container. In view of the 
comments made above, we see that formulation of an optimal control problem for the 
original time-dependent system would require us to satisfy certain optimality conditions 
on the boundary of each individual moving droplet in addition to conditions on the free 
surface of the liquid in the pool. On the other hand, the concept of the effective free 
surface allows us to replace this optimization problem with a simpler one, which is also 
computationally more tractable, where the optimality conditions have to be imposed on 
the stationary effective surfaces. Thus, as one application, the proposed approach will 
allow us to extend the optimization formulation developed in Volkov et al? to include 
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FIG. 1. Schematic of our model problem representing droplets impinging on a free surface, a 
phenomenon typically encountered in various welding processes such as Gas-Metal Arc Welding 
(GMAW). The solid lines represent the actual time-dependent interface between the liquid and 
gas phases, whereas the dashed line is the steady effective surface Ylg we seek to determine. 



the effect of the mass transfer into the weld pool via droplet impingement. 

We remark that droplet impingement onto a thin liquid film is a phenomenon with 
manifold manifestations in technology, including chocolate processing, spray painting, 
corrosion of turbine blades, fuel injection in internal combustion engines, and aircraft 
icing. It also occurs in many natural phenomena such as the erosion of soil and the 
dispersal of spores and micro-organisms. A considerable amount of literature is available 
as concerns the numerical modeling of droplet impingement onto a solid surface. Harlow 
& Shannon 28 were the first to simulate this phenomenon and several other authors have 
applied the Volume of-Fluid (VoF) based approaches such as RIPPLE 29 and SOLA- 
VOF 30 to understand droplet impingement phenomena. Trujillo et al. 31 also performed 
a numerical investigation and experimental characterization of the heat transfer from a 
periodic impingement of droplets. 

The structure of this paper is as follows: in the next Section we present the formu- 
lation of the problem in general terms, in the following Section we introduce our model 
problem and in Section IV we discuss a very simple closure strategy which may be suit- 
able for this problem, in Section V we introduce a shape-optimization approach to the 
numerical solution of the resulting averaged equations, whereas in Section VI we present 



4 















V / [ I 

I \ 11 

* \r 


if)/ G 


S lg 

\ \ \ 


* \ x 

\ V 
\ \ s 




'Jl 



FIG. 2. Schematic of the domains and domain boundaries used in the definition of the model 
problem in Section II A. The domain Vl^if) occupied by the liquid phase is marked in gray, 
whereas the thin and thick solid lines represent, respectively, its boundary Y^cif) and the 
corresponding effective surface Tlg- The subregion Vt r (see Section II C) is delimited by the 
thick dashed lines. 

some computational results together with a discussion; final conclusions are deferred to 
Section VII and some technical results concerning solution of the shape optimization 
problem in Section V are collected in Appendix A. 

II. PROBLEM FORMULATION 

In order to simplify the presentation of our approach, we will consider a two- 
dimensional problem formulated in a general domain Vt C K 2 , shown schematically 
in Figure 2, where Ql and Qq represent the subdomains occupied, respectively, by the 
immiscible liquid and gas phases, whereas Y L g represents the liquid-gas interface (e.g, 
droplet boundary or the free surface of the weld pool). 

A. Assumed Governing Equations 

For a general description of the equations and boundary conditions governing mul- 
tiphase flows we refer the reader to monograph by Prosperetti & Tryggvason 27 . We 
assume that our model problem involves the following dependent variables 

(a) velocity v = [u, v] T : Q x (0, T] M 2 , 
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(b) pressure p : fix(0,T] -> R and 

(c) position of the free surface V tG (o,r], ^LG(t) — ^l(^) H ficW, 

where T denotes the length of the time window of interest and "=" means "equal to 
by definition". It is also assumed that there is no mass transfer across the interface 
Tlg- We have the following equations governing the fluid flow in the two phases, for the 
moment in the time-dependent form 

PL-^ + Pl(v • V)v - V • a - p L g = in Sl L , (la) 

V-v = inft L , (lb) 



<9v 

pG ~dt + PG<(V ' ~ V a ~ PGg = ° in ^ 

V • v = in tt Gl (2b) 

where Pl and pc are the densities in the liquid and gas phase and a = —pi + cr^ is the 
stress tensor in which cr^ = /x[Vv + (Vv) T ] , I denotes the identity matrix and the vis- 
cosity coefficient p = Pl or p = pc m the liquid and gas phase, respectively. The symbol 
g denotes the gravitational acceleration. Equations (lb) and (2b) represent conservation 
of mass, whereas equations (la) and (2a) represent conservation of momentum in both 
the liquid and gas phase. Systems (1) and (2) are subject to the following boundary 
conditions on the liquid-gas interface T L g 

v| L = v| G on T LG , (3a) 

G 

.n = 7^n on Y LG) (3b) 

L 

where n and t are the unit normal and tangential vectors on the interface T L g-> k is 
the interface curvature, 7 the surface tension (a material property assumed constant), 
whereas the subscripts L and G (with or without the vertical bar) denote quantities 
defined in the corresponding phases (Figure 2). We note that the vector-valued condition 
(3b) implies the balance of both the normal and tangential stresses. For simplicity, on 
the far boundary r , cf. Figure 2, we adopt the no-slip boundary condition 

v| G = on r . (4) 

As regards the mathematical description of free-boundary problems, there are two 
main paradigms, namely, (i) "interface tracking" approaches, see Neittaanmaki et al 7 
and (ii) "interface capturing" approaches, see Sethian 8 . While description (l)-(3), fea- 
turing the location of the interface Y L g as the dependent variable, belongs to the first 
category, for the purpose of developing our formulation an interface capturing approach 



6 



will be more suitable and we employ a technique known as "Volume of Fluid" (VoF). 
However, our computations of the effective surfaces will be ultimately carried out with an 
"interface tracking" approach, see Section V. A detailed description of the VoF method- 
ology can be found in the paper by Hirt & Nichols 30 , see also monograph by Prosperetti 
& Tryggvason 27 . This method employs the "volume fraction" as an indicator function 
to mark different fluids 

While in the continuous setting the interface Y^g is sharp and the VoF function F may 
take the values of and 1 only, in a numerical approximation there may exist a transition 
region where < F < 1 and the fluid can be treated as a mixture of the two fluids on 
each side of the interface. The values of the indicator function are associated with 
each fluid and hence are propagated as Lagrangian invariants. Therefore, the indicator 
function obeys a transport equation of the form 

dF 

— + (v • V)F = in fl (6) 

dt 

Based on the indicator function, local material properties such as the density p of the 
fluid can be expressed as 

p(F(x)) = F(x)p L + [1 - F(x)]pg. (7) 

Relationship (7) allow us to rewrite formulation (l)-(3) in an equivalent form as one 
system of conservation equations defined in the entire flow domain Q where the fluid 
properties are, in general, discontinuous across the interface between the two fluids. In 
this single-field representation the two fluids are identified by indicator function (5), 
whereas the material properties are expressed as piecewise constant functions and can 
be written in terms of their values on either side of interface T LG: cf. (7). 

dp ^ + V • [p(F) v] = in 0, (8a) 



dt 
dp(F) v 

dt 



+ V • [p(F) vv] = V • a + / n<5(x - x') ds(x') in 0, 

•'I;..,- 



(8b) 



d J- + ( v • V)F = in Q, (8c) 

where vv denotes the dyadic product, i.e., the tensor defined as [vv]^. = [v^ [v] 
z,j = 1,2. Conservation equations (8a) and (8b) can be obtained in a straightfor- 
ward manner by considering the integral balance of mass and momentum for the fluid 
with variable density p(F) in some arbitrary control volume. Further discussion of the 
single-field description of two-phase flows can be found in the monograph by Prosperetti 
& Tryggvason 27 . 
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The last term on the right hand side (RHS) in (8b) represents the source of momen- 
tum due to the surface tension. It is related to boundary condition (3b) and only acts 
at the interface Y L g as indicated by the presence of the Dirac delta function in the inte- 
grand expression of the integral. The surface integral in equation (8b) can be difficult to 
evaluate directly. In order to overcome this problem, a continuum surface force (CSF) 
model was introduced by Brackbill et al 32 which represents the surface tension effects in 
terms of a continuous volumetric force acting within the transition region which arises 
when the problem is discretized. The surface integral in (8b) is therefore approximated 
as in Prosperetti and Tryggvason 27 

/ 7*/ n'5(x - x')dT « VF, (9) 

whereas the curvature of the interface can be computed in terms of the VoF function as 
follows 



K=v 'lr^J- (10) 

Using (9) in (8b), we obtain a simpler form of the one-field system (8), namely 

dp(F) 



dt 

dp(F) v 



V-[p(F)v]=0 in ft, (11a) 

+ V • \p(F) vv] = V • a + 7^ VF in ft, (lib) 



dt 

dF 

— + (v • V)F = in ft. (11c) 

dt 



B. Averaging Procedures 

The goal of this Section is to derive a time averaged form of governing system (1)- 
(3), or equivalently (11), and state the "closure problem", i.e., identify the terms in the 
resulting equations which need to be modeled. Our objective is to express the averaged 
equations solely in terms of averaged velocity, averaged pressure and averaged indicator 
function as the dependent variables. A number of different averaging techniques have 
been considered in the literature in regard to multiphase flows 17 ' 18 ' 20 ' 25 ' 27 . Here we will 
rely on the conventional time-averaging procedure, see Monin & Yaglom 9 which is based 
on the ideas originally due to Reynolds (it should be added that in statistical physics 
averaging is typically performed with respect to realizations, however, in view of the 
ergodicity assumption adopted here, the ensemble average can be replaced with a time 
average used in (12) below). Given the quantity : [0, T] x ft —> ]R d , d — 1, 2, we thus 
define the pointwise time average as 

i rto+At 

M( x ) 4 At/ t ^' x )^' ( 12 ) 
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where the time window At is assumed large compared to the time scale of the random 
fluctuations associated with free boundaries. Since in the present problem we are in- 
terested in steady solutions, we take At —> oo, so that the averaged variables do not 
depend on time. In the conventional Reynolds decomposition, the chaotically varying 
flow variables are replaced by the sums of their time averages and fluctuations, i.e., 

v = (v)+v', p=(p)+p', p =(p)+p>. (13) 

By definition, the time average of a fluctuating quantity is zero, i.e., (V) = 0, (//) = 
and (p f ) = 0. We also note that the averaging operator (•) commutes with differentiation 
with respect to the space variables 9 ' 17 . We shall furthermore assume that 9 

Our derivation of the averaged equation follows the general development presented 
in Hong & Walker 18 , although we use a somewhat different notation adapted to the 
present problem. We begin with continuity equation (8a) and decompose the dependent 
variables as in (13). The equation is then time-averaged and we obtain 

V-«P><v» = -V-« P V)). (15) 

We need to re-express the right hand side of equation (15) to eliminate p' . From (7) 
and applying the Reynolds decomposition to the indicator function (5) 

F(t,x) = (F)(x) + F'(t,x), where (F')(x) = 0, (16) 

we obtain 



(17) 



p(t,x) = [(F) (x) + F'(t,x)]p L + [1 - <F)(x) - F'(t,x)]p G 
= (p) + F'(t,Tt)(p L -p G ) 
which allows us to identify p'(^ x ) = F'(t,x) (p L — pc). Using (17) we can now deduce 

(pV) = (p L -p G )(FV), (18) 

so that (15) becomes 



V 



{ [ PL (F)(s) + PG (1 - <F)(x))] <v»} = -(p L - p G )V • <FV>. (19) 



which is the Reynolds-averaged form of the continuity equation, where the right-hand 
side (RHS) terms represent the average effect of the fluctuations of the free boundary. 

We now turn our attention towards momentum equation (11). In order to simplify 
the formulation of the present problem we make the following 

Assumption 1. The fluctuations of viscosity /x(i, x) = // L F(t, x) + /xg[1 — F(£, x)] and 
the interface curvature cf (10) ; are neglected. These quantities will be therefore treated 
as constant and will not be subject to averaging. 
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The reason for this simplification is that proper handling of viscosity and curvature 
fluctuations leads to significant complications in the resulting average equations, and 
this issue is deferred to future research. In a phenomenon characterized by an interplay 
of capillary, viscous and inert ial effects, Assumption 1 implies that the applicability of 
the model is restricted to flow regimes dominated by the inertial effects. Indeed, we 
expect that the density fluctuations are going to have a dominating influence on the 
effective surfaces in the class of applications motivating this study. The development of 
the Reynolds-averaged form of the momentum equations proceeds most easily when the 
nonlinear advection terms are written in the conservative form. Again, at every point x 
we replace the dependent variables with relations (13) and then average the equations 
over time. The complete Reynolds-averaged momentum equation can be written as 



V-((p)(v)(v)) = -V( P )+V-(K)-(p)(vV)-(v)( P V)-(p'vV))+7kV(F), (20) 



where (a* 1 ) is the usual viscous stress tensor defined in terms of the averaged velocity 
field and, in addition to the Reynolds stresses, on the RHS in (20) we also note the 
presence of new terms representing fluctuations of the free boundaries. 

The main idea behind the proposed approach is that the resulting averaged solutions 
should preserve some essential features of the original time-dependent free-boundary 
problem (l)-(3), namely, a sharp separation between the two phases, cf. Figure 1, along 
an interface which we defined as the effective free surface. While the time-dependent 
indicator function F may only assume the values of and 1, cf. (5), its average (F) may 
assume all intermediate values < (F)(x) < 1 (this is in fact clearly visible in the plots 
of the mean indicator function (F) obtained by averaging the solutions of our model 
problem, see Figure 4a to be discussed further below in Section IV). Such smoothly 
varying indicator functions (F) correspond to a continuous transition between the two 
phases without a well-defined interface. Therefore, in order to be able to define averaged 
flows with sharp effective boundaries we have to introduce the following 

Assumption 2. The average indicator function (F) is replaced in Reynolds-averaged 
equations (19)-(20) with the piecewise-const ant function F : Q — >• K. defined as follows 




(21) 



where and Qg are the corresponding time-invariant subdomains occupied by the liquid 
and gas phases. 



With this assumption the Reynolds-averaged equations take the form 
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where the vector A and tensor B are defined as 



Pg) 
Pg) 



(F'u') 
(F'v') 

2{u){F'v') + (FW) (u)(F'v') + {v){F'u') + (FW) 
(u) (F'v f ) + (v) {F'u') + (FW) 2(v) (F'u') + (FW) 



(23a) 



(23b) 



As one can see, equations (22) are not "closed" , because they contain averaged products 
of fluctuation terms for which no additional equations are available. Therefore, we 
will seek to model these terms with closure expression of the form A = A(F, (v)) 
and B = B(F, (v)) which are functions of the averaged dependent variables F and 
(v) = [(u), (v)] T '. This is in addition to the closures required for the "classical" Reynolds 
stress tensor (vV). While modeling the latter expressions is a rather well-advanced 
area 6 , development of closures for product terms corresponding to fluctuations of the 
free boundaries has been the subject of relatively few investigations, see, e.g., Refs. 20 
and 25, which focused on the case with a diffuse interface. A very simple closure model 
for these terms adapted to the present formulation of the problem with a sharp interface, 
cf. Assumption 2, will be presented in Section IV. The question of closure models for 
the Reynolds stress terms will not be considered in this work. 



C. Reduction of Averaged Fluctuation Terms to Boundary Conditions 

In the derivation of the closure models the quadratic and cubic products involving 
the fluctuation fields F', u' and v' will need to be expressed solely in terms of the 
time-averaged fields F, (u) and (v). As regards the dependence on F, this means that 
expressions for these closures will depend on the location relative to the effective free 
surface and, evidently, the components of the tensors A and B are nonvanishing only 
in a close proximity of the free boundary TlGi cf. Figure 2. From the point of view 
of the formulation of a computation-oriented model it is therefore not "economical" 
to introduce new terms into the averaged equations which would be nonzero only in 
a very small fraction of the domain. We therefore propose the following simplifying 
approach in which the averaged fluctuation terms involving tensors A and B defined in 
the bulk are approximated with suitable terms defined on the effective boundary Tlg- 
This can be done by integrating the terms involving A and B in (22a)-(22b) over their 
support Q f = supp A = suppB and then using the divergence theorem (in principle, the 
supports of these two terms may in general be different, but for the sake of simplicity 
we assume here that they coincide; this simplification does not in any way affect the 
accuracy of the proposed approach) . We remark that analogous ideas were also pursued 
by Brocchini & Peregrine 20 and by Brocchini 21 . One important difference between these 
approaches and the formulation explored here concerns the description of the effective 
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boundary (explicit in Refs 20 and 21 versus intrinsic considered here). Noting that the 
fields A and B are discontinuous at the effective surface Y L g (which is contained inside 
the integrations domain Q J ') ) and vanish on dVt 1 we obtain 

h = I V-AdQ= I [xi'A] G L da=[ a da, (24a) 
Jw Jf LG Jr LG 

I 2 = [ V-Bdtt=[ [n-B]^ da=[ bda (24b) 
Jw Jf LG ° Jr LG 

in which the fields a : f L q —> K and b : f L q —> K 2 are defined in terms of the jumps 
of A and B as 



= [n-B]?. (25) 



We thus see that in the mean sense the fluxes due to the fluctuating terms V • A and 
V • B in the averaged mass and momentum equations (22a) and (22b) can be realized 
by the terms a and b, cf. (25), defined on the effective boundary Tlg- This leads to the 
following 

Assumption 3. which has two parts 

(a) we replace the source term V • A in averaged mass conservation equation (22a) 
with an additional term (an) in the corresponding boundary condition (3a), 

(b) we replace the source term V • B in averaged momentum conservation equation 
(22b) with an additional term b in the corresponding boundary condition (3b) , 

so that the following system of equations is obtained (rewritten here in the two subdo- 
mains together with all boundary conditions) 



V) (v) - V • 


(cr) - p L b 


= 


in Ql, 


(26a) 




V-(v) 


= 


in fli, 


(26b) 


V)(v) - V- 


(cr) - p G g 


= 


in Qq, 


(26c) 




V-(v) 


= 


in Q,q, 


(26d) 




[<v)£ 


= an 


on r LG , 


(26e) 




M(cr)g 


= 7K n + b 


on r LG , 


(26f) 


«v>L + 


<v)| G )-n 


= 


on f LG , 


(26g) 



where boundary condition (26g) corresponds to condition (4) in the situation when the 
normal velocity at the effective surface is allowed to have a discontinuity, cf (26e). 

As is evident from Figure 2, this Assumption is satisfied when the subregion Q' forms 
narrow bands along the effective free boundary Y L g which happens when the fluctua- 
tions of the free boundary T L g occur at a length-scale significantly smaller than the 
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characteristic dimension of the entire domain. As regards the averaged conservation 
equations, Assumption 3 has the effect of reducing, or localizing, the influence of the 
averaged terms involving fluctuation of the free boundary to the effective free boundary 
f L q. The as of now undefined functions a = a(F, (v)) and b = b(F, (v)) represent the 
required closure models and need to be determined separately for every flow problem. 
We add, that since these functions depend on the location of the effective free surface 
Y L g ) boundary conditions (26e) and (26f) are in fact geometrically nonlinear. We also 
remark that in Ref. 20 closure models for certain free boundary problems were derived 
based on an analogous concept of integral balances in the surface layer. Construction 
of a very simple closure model for functions a and b applicable to a model problem 
introduced in the next Section will be presented in Section IV. 

III. MODEL PROBLEM 

While up to this point our discussion has been concerned with a generic two-phase, 
free-boundary problem, we will from now on focus on a specific flow configuration. 
Thus, to fix attention, we will consider the flow set-up shown schematically in Figure 
1. It features droplets entering the domain periodically through the top boundary 
and impinging on the free surface resulting in sloshing. On the lateral boundaries T 
no-slip boundary condition (4) is applied and we observe that, respectively, an unsteady 
or steady contact line will appear where the time-dependent interface T^, or the corre- 
sponding effective surface Y L g ) intersects the boundary r . While it is well known that 
subject to the classical no-slip and free-surface boundary conditions the contact-line 
problem is not well-posed 33 , development of a both mathematically and physically con- 
sistent description of this problem still remains an open question. Addressing this issue 
is beyond the scope of the present investigation, and our treatment of the contact line is 
a standard one: in the solution of the time-dependent problem a suitable regularizing 
effect is achieved by discretization of the governing equations (described further below), 
whereas in the solution of the steady problem with the effective surface regularization 
is introduced via formulation in terms of variational shape optimization. Application 
of this numerical approach to a closely related problem with a contact line singularity 
is analyzed in detail by Volkov and Protas 10 . In order to maintain a constant average 
(over one period of droplet impingement) mass of the fluid M = J Q F(t, x) dQ, the fluid 
is drained through the bottom boundary of the domain (i.e., suitable nonzero velocity 
boundary condition v • n 7^ is applied there) . 

Solutions to the problem described above depend on the following three parameters 

(a) length T of the interval at which droplets are released, 

(b) velocity of the droplet V^, and 

(c) radius of the droplets r. 
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We emphasize that the choice of this particular model problem was in fact inspired by an 
industrial application described in Volkov et. al? which has also motivated our broader 
research program. Numerical solution of this time-dependent free-boundary problem 
is obtained using the solver InterFOAM which is a part of the library 0penF0AM 13 and is 
based on the VoF method. Details of the numerical method and its implementation in 
InterFOAM can be found for instance in the Ph.D. thesis of Rusche 14 . The resolution 
used in our calculations was 100 x 100 grid points with a nondimensional time step of 
0.05. In order to characterize the time-dependent and mean fields obtained as solutions 
to this problem, in Figure 3 we present several snapshots of the indicator function F(£, x) 
at different time levels spanning two periods of droplet impingement. To fix attention, 
the results presented in Figure 3 were obtained using the following parameters T = 1.0, 
V d = 1.0 and r = 0.25. 

IV. ALGEBRAIC CLOSURE MODEL 

In Section II B we introduced the Reynolds decomposition of the flow variables into 
the time-averaged quantities (denoted with angle brackets (•)) and fluctuating quantities 
(denoted with primes). The terms involving averaged products of fluctuating quantities 
appear as unknowns in averaged equations (22) and must be closed with suitable "closure 
models", analogous to those which arise in classical turbulence modeling approaches. 
The most commonly used methods of turbulence modeling are surveyed in monograph 
by Pope 6 . Briefly speaking, depending on their mathematical structure, such approaches 
fall into two main categories, namely, algebraic models and differential models in which 
evolution of the quantities introduced to close the system is governed by additional 
PDEs. Some attempts at deriving closure models for two-phase flows were already made 
by Brocchini & Peregrine 20 and by Brocchini 21 who obtained such models for regimes 
characterized by different values of the turbulent kinetic energy and the turbulent length 
scale. In this Section, we make an attempt to derive an extremely simple algebraic closure 
relationship based on an elementary model of the process defined by the following set of 
assumptions (see also Figure 4). 

Assumption 4. (a) Droplets are spherical with radius r and move as rigid objects, 

(b) there is no collision or coalescence of droplets, 

(c) droplets are falling periodically with frequency T _1 and constant velocity Vd, 

(d) the fluid outside droplets (i.e., the gas phase) is motionless, 

(e) the mean fields do not depend on the vertical coordinate. 

We observe that Assumption 4b constrains the problem parameters so that 2r < Vd T. 
It is also to be noted that Assumption 4e implies that the model is effectively one- 
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(g) t = (3/2)T (h) t = (7/4)T (i) * = 2T 

FIG. 3. Snapshots of the indicator function x) obtained at the indicated instants of time 
in the solution of the time-dependent problem (l)-(4) with the parameters T = 1.0, = 1.0 
and r = 0.25. 

dimensional with variations only in the direction normal to the effective surface. While 
the above assumptions are rather far-reaching (in particular, the model does not include 
any effects of droplet impingement on the free surface), our objective here is to provide 
some preliminary insights concerning computation of effective free surfaces, and develop- 
ment of closures based on more accurate models is left to future research (some possible 
directions are discussed briefly in Section VII). We thus proceed to use Assumptions 4 
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(a) (b) 

FIG. 4. Sketch illustrating the main features of the model based on Assumptions 4: (a) 
construction of the piecewise-constant approximation F to (F), (b) time-dependence of the 
indicator function F(x,t) for different values of the distance x. 

in order to derive expressions for the fluctuating fields F' ) v! and v' which will be given 
in terms of the mean fields (F) (or F), (u) and (v). These expressions will be in turn 
used to determine the fields a and b in (26e)~(26f). 

The coordinate system is shown in Figure 4. We begin by observing that in the 
model problem considered the horizontal velocity component vanishes identically, i.e., 
u(t, x) = 0, t > 0, x G K, as do its mean and the corresponding fluctuation fields 

(u)(x) = 0, u'(t, x) = 0, t > 0, x G R. (27) 

Since the model considered assumes periodic behavior, without loss of generality we 
are going to focus on a single period of droplet impingement, i.e., t G [0,T], and we 
also remark that the vertical velocity component v and the indicator function F are 
piecewise-constant functions of time at every point in space. We thus define the following 
"pulse" function 

^ ^ o, otherwise ' ^ ^ 

where < 9 < T, which allows us to write the following expression for the indicator 
function F as a function of time and the coordinate x (for simplicity, we omit the y- 
dependent phase shift in this expression, as it does not affect the averages which we will 
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ultimately compute) 

f n 2vCT (t), < |x| < r 

F(t,x)=\ ' , (29) 

[0, \x\>r 

and then the vertical velocity component becomes v(t, x) = Va F(t, x) for t > 0, x G R. 
Next, computing the average over one period we obtain 

j, ( 2Vr 2 - x 2 

(F)(x) = ±f F(t,x)dt=l V d T ' ^\ x \<\ (30) 

I 0, |x| > r 

and also have (v)(x) = Va (F)(x). These expressions allow us to evaluate the fluctuating 
fields as follows 



2\/r' — 



V d H 2 _v^(t)-^; , < |^| < r 



" 2 rjQ 2 

v'(t,x)=v(t,x)-{v)(x) ={ —2^^ T ' ^, (31) 

0, |x| > r 



f (t, X ) = F(i, x) - (F)(x)= { "^F ( " 'vj 1 ' -' 11 ''. (32) 

0, \x\ > r 



We note that for — r < x < r the averaged indicator function (F)(x) has a quadratic 
distribution which can be interpreted as resulting in a smeared interface. However, in 
view of Assumption 3, we require a sharp interface Tlg corresponding to a piecewise- 
constant indicator function F 





X 




X 



F = < „ > (33) 



where r is the new interface location which can be determined based on the principle 
of mass conservation. It is expressed using the original smeared (30) and the new 
piecewise-constant (33) indicator functions as follows 

POO POO 

/ (F)(x)dx= / Fdx (34) 
Jo Jo 



from which we obtain 



7rr 2 



ro = 2^T- (35) 
In view of Assumption 4b, it is evident that Tq < r. Therefore, one can recalculate the 
fluctuating indicator function, now with respect to F given by (33) and (35), to obtain 

HsxAaW - 1 I x l< r o, 



F\t ) x)={ n 2AAa (t) r <|x|<r, (36) 

I x |> r. 
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We remark that the obtained expressions (31) and (32) for the fluctuating quantities 
depend only on the model parameters {T, r, Vd} and the position of the effective (sharp) 
interface t L c- We are thus in the position to calculate the averaged products of fluctua- 
tions appearing as components of tensors A and B in (23a) and (23b). First, we observe 
that in view of (27) we have 

(FV) = (F'u'u') = (F'u'v') = (37) 

everywhere in the domain. As regards the products of fluctuations which do not include 
u', we observe that the form of the expression will depend on the coordinate x, and the 
following three regions are distinguished 

(a) inner region defined by | x |< r , i.e., where F(x) = 1 and (F)(x) > 0, 

(b) transitional region defined by r <| x |< r, i.e., where F{x) = and (F)(x) > 0, 

(c) outer region defined by | x |> r, i.e., where F(x) = (F)(x) = 0. 

The expressions for the averaged products of fluctuating quantities in these regions are 
discussed in three subsections below, whereas in the last subsection we demonstrate 
how these expressions can be used to close the terms a and b in boundary conditions 
(26e)-(26f). 



A. Inner Region 

In the inner region \x\ < ro, using (31) and (32), we can write after employing some 
straightforward properties of pulse function (28) 



Time-averaging expression (38) we obtain 



(Fv) = — (39) 



Following the same steps as above, we can deduce that 



r 



- x 2 ) 3 / 2 4(r 2 - x 2 ) 



18 



B. Transitional Region 



In the transitional region r < \x\ < r, again using (31) and (32), and following the 
steps involved in obtaining (39) and (40), we can express (FV) and (F'v'v') as 



(F'v') = 



2\/r 2 - x 2 A(x 2 - r 2 ) 

t vyr 2 ' 



2V d Vr 2 - x 2 8(r 2 -x 2 ) 8(r 2 - x 2 f/ 2 

(Fw) = Vd T* ■ 

C. Outer Region 

Since in the outer region \x\ > r we have F'{x) = 0, it follows that 

(F'v f ) = 0, 
(F'v'v') = 0. 



(41a) 
(41b) 



(42a) 
(42b) 



D. Closure Terms in the Boundary Conditions on the Effective Surface 

In this Section we derive expressions for the closure terms a and b in boundary 
conditions (26e)-(26f) based on the simple algebraic closure model proposed here. First, 
since u' = and (u) =0, we note that 



A = (p L - Pg) 





(F'v') 



B = (p L - p G ) 





(F'v'v') 



(43) 



where the form of the nonzero entries depends on the location with respect to the 
effective surface T L g (see Sections IV A, IV B and IV C). In our model the location of 
the effective surface Y^g coincides with the boundary between the inner and transitional 
regions, cf. Figure 4a. In view of (25) and (43) we therefore have 



a = (p L -p G ) {n y {F'v')\ 
b = (pl ~ Pg) 

= (Pl ~ Pg) 



transitional 



n 















n y (F'v'v') 




transitional 


n y (F'v'v') 





T 2 



^4V 2 T 2 -7v 2 r 2 - ^ (4V 2 T 2 - n 2 r 2 ) 



(44) 



(45) 



where we also used the relationships derived in Sections IV A and IV B evaluated at 
x = r to express the fluctuating terms. Relations (44) and (45) close our system of 
averaged equations (26). The conclusion that a = is consistent with the fact that 
there is no mass production at the effective surface. It should be emphasized that 
the proposed closure model is quite problem-specific and it is not obvious whether the 
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same ideas could be used to develop closures for other flows with effective free surfaces. 
The closure model derived by Brocchini & Peregrine 20 for an analogous flow regime is 
more general, and at the same time less explicit, as it is constructed in terms of an a 
priori unspecified "probability function" describing ejection of droplets with some specific 
velocity. The model proposed in this Section can be regarded as assuming a specific form 
of this probability function. A numerical approach well suited for the solution of free- 
boundary problem (26) is described in the next Section, whereas in Section VI we will 
analyze the dependence of the solutions on the parameters {V^, i?, T} characterizing the 
closure model. 

V. SOLUTION OF AVERAGED EQUATIONS WITH EFFECTIVE 
SURFACES VIA A SHAPE OPTIMIZATION APPROACH 

In Sections II and IV we formulated a set of steady-state PDEs as a simplified 
time-averaged model of a fluid problem involving unsteady free boundaries such as, 
for example, the system introduced in Section III with droplets impinging on the pool 
surface, cf. Figure 1. In the proposed formulation, the steady liquid-gas interface is 
represented as the effective free surface described mathematically as the discontinuity of 
the averaged indicator function F. Therefore, the set of governing equations (26) has the 
form of a steady free-boundary problem. Since such problems tend to be hard to solve 
numerically, we argue below that a computationally efficient approach can be developed 
by formulating this problem in terms of shape optimization. More specifically, we will 
frame it as finding an optimal shape of the interface (i.e., effective free boundary T L g) 
such that a cost functional representing the residual of one of the interface boundary 
conditions will be minimized with respect to the position of the interface subject to 
the constraints representing the governing (time-averaged) equations and the remaining 
boundary conditions. We refer the reader to the monograph by Neittaanmaki et al 7 
for a general discussion of advantages of such an approach, and to papers by Volkov & 
Protas 10 and Volkov et al 2 for a discussion of some applications to problems similar to 
the one studied here which also included treatment of contact lines. 

System (26) represents a steady-state free-boundary problem where T LG has to be 
found as a part of the solution. By fixing the domain and its boundary Tlq, and 
removing one of the boundary conditions, for example, the normal component of (26f), 
we obtain a steady fixed-boundary problem. The residual of the normal component of 
condition (26f) is then minimized with respect to the unknown shape of the effective 
surface T L g using a suitable shape-optimization algorithm. The choice of the normal 
component of condition (26f) as the optimization criterion is motivated by the fact that 
this condition contains closure terms, hence in practical situations need not be satisfied 
up to the machine accuracy. 

Our solution approach is then formulated as follows. Suppose we define the function 
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FIG. 5. Schematic of a single step of the shape optimization algorithm (50): (solid line) 
approximation to the effective surface at the fc-th iteration and (dashed line) approximation 
to the effective surface at the k + 1-th iteration. 



X : ^lg R which will represent the residual of the normal component of the 
momentum boundary condition (26f) 



X = n • (<r£) • n - n • (<r£) -n-jK - b • n. 
The cost functional can then be defined as 

J$lg) = \( X 2 da, 
so that the optimization problem becomes 

mmJ{t LG ), 



(46) 



(47) 



(48) 



where x m (46) depends on the shape of the effective free surface which is the 

control variable in optimization problem (48), via governing PDEs (26a)~(26d) subject 
to boundary conditions (26e), (26g), and the tangential component of condition (26f), 
i.e., 



n 



0, 



on 



LG- 



(49) 



Thus, we observe that the system of PDEs serving as the constraint for optimization 
problem (48) has in fact the form of a fixed-boundary problem which makes evaluation 
of the cost functional at every iteration easier (i.e., it corresponds to the approximation 
T^q of the effective surface at the given k-th iteration). The position of the effective 
interface Y L g can then be found using the following iterative gradient-descent algorithm 



xL(fc+i) = xL(fc) +r k G Vjit^ 

1 LG 1 LG L V 



LG 



n 



k 



1,2, 



(50) 
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ALGORITHM 1: Iterative minimization algorithm for solving system (26) via a shape 

optimization approach. 
Input: e r (adjustable tolerance), (initial guess for the effective surface) 
Output: Tlg (effective surface) 

n <- 

T^Q ^— initial guess (Figure 6) 
repeat 

solve direct system 

solve adjoint system 

compute gradient VJ (v 

perform line minimization to determine 

update effective surface using (50) 

n <— n + 1; 

until rW < e r . 



where xL(fe) G Mr represents points on the interface Tlq at the fc-th iteration and t& is the 

1 LG 

length of the step in the descent direction. The function G determines the specific form 
of the optimization algorithm used (e.g., the steepest descent, conjugate gradients, or 
quasi-Newton method, etc., see Ref. 15). In our results reported in Section VI we use the 
Conjugate Gradients Method. A central element of algorithm (50) is the cost functional 
gradient V 1 J : T y L Q — >• M representing the continuous sensitivity of cost functional 
(47) to infinitesimal modifications in the normal direction of the shape of f^. In 
other words, as indicated in Figure 5, the scalar-valued function V J [v , depending 
on the arclength coordinate along the interface, represents the normal displacement of 
the current approximation to the effective surface f ^ resulting in the largest possible 
decrease of cost functional (47), see also Appendix A 2. The contact points, where 
the effective surface T L g meets the solid boundary r , require special attention and we 
follow here the approach developed by Volkov and Protas 10 to deal with the contact 
line problems in the context of shape optimization. Determination of the gradient V J 
requires the solution of a suitably-defined adjoint system and details concerning its 
derivation are deferred to Appendix A. The step size t& is obtained via solution of the 
following line minimization problem 

r k = argmin r>0 J (x| f ^ + r G [v J (f g}) n] ) (51) 

which can be done, for example, using Brent's method, see Ref. 15. The complete 
approach is summarized in Algorithm 1. 

Our time-dependent model problem (l)-(4) is set up such that the mass of liquid 
remains constant over every period of droplet impingement (the amount of mass drained 
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at the bottom of the container equals the mass added in the form of droplets at the top of 
the domain, cf. Section III). Therefore, it is necessary to ensure that in our steady-state 
averaged problem with effective surfaces (26) the same mass is enclosed in the liquid 
domain Q^. Mathematically, this is implemented by constructing an initial guess for 
the liquid domain which has the prescribed mass and then making sure that this 
mass is not changed at subsequent iterations. For this to happen, it is required that the 
shape gradients VJ{T L q) do not change the volume of Ql. This property is enforced 
at each iteration as follows. First, we calculate the mean value M G K of the gradient 
on the effective surface 

M = | / VJ(s)ds (52) 
L Jo 

where s G [0, L] is the corresponding arclength coordinate. The new gradient with zero 
mean displacement in the normal direction is then obtained as 

VJ(s) = VJ(s)-M, V sG[0 ,l]. (53) 

The cost functional gradient V J is replaced with zero-mean gradient V J in expressions 
(50) and (51). 

VI. RESULTS AND DISCUSSIONS 

In this Section we present sample computations for the problem of determining the 
effective free surfaces in the flow described in Section III, see also Figure 1. In order 
to calculate the cost functional gradient (given by expression (A18) in Appendix A) 
we need to solve "direct" system (26) and adjoint system (A13)-(A14). Both these 
solutions are obtained using the finite-element method implemented in the COMSOL 
script environment. The domain (Figure 6) is discretized using approximately 4000 
Lagrangian elements with mesh size varying between 0.9 to 0.01. In all computations 
presented here we used the physical parameters with values indicated in Table I and 
these calculations were performed using the Navier-Stokes and Poisson solvers available 
in COMSOL. For illustration purposes, in Figures 7a and 7b we show the fields of the 
direct and adjoint vorticity obtained at the first iteration. Before analyzing the solutions 
with effective free surfaces obtained for different parameters of the closure models, we 
validate the calculation of the cost functional gradient which is the main element of our 
computational approach, cf. Section V and Appendix A. 

A. Validation of the Shape Optimization Approach 

In this Section we demonstrate the consistency of the gradients ^y lg J obtained 
using expression (A18) in Appendix A3. A standard test consists in computing the 
Gateaux differential of cost functional J(T L g) in some arbitrary direction x' = ('n 
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Physical Parameter 


Value 


Density of the liquid, pl 
Dynamic viscosity of the liquid, pl 
Density of the gas, pc 
Dynamic viscosity of the gas, pc 
Gravitational acceleration, g 
Surface tension, 7 


1000[i^.m- 3 ] 
0.001[Pa.s] 
l[Kg.m~ 3 } 

0.00001 [Pa.s] 
9.81[m.5- 2 ] 
O^.m' 1 ] 



TABLE I. Values of the physical parameters used in the computation. 




-2 5 -2 -1.5 -1 -0.5 



FIG. 6. Geometry of the computational domain £1 with the liquid and gas subdomains (l^ and 
Qg and the effective boundary Tlg used as the initial guess in shape optimization algorithm 
(50). The figure also shows the unstructured triangular finite-element mesh used in the solution 
of direct and adjoint problems. 

using a finite difference technique and comparing it with the expressions for the same 
differential obtained using the gradient ^y lg J and Riesz representation formula (A7). 
The ratio of these two expressions, which is a function of the finite-difference step size 
e, is defined as 

J(x f + <-n)-J(f LC ) 

Proximity of ^f LG ( 6 ) ^° ^ ne uni ty is thus a measure of the accuracy of the cost functional 
gradient computed based on the adjoint field. Figure 8 shows the behavior of the quantity 
function of the parameter e for different perturbations (' . We note that in 
all cases the quantity ^f LG ( e ) * s Q ui te close to unity for e spanning over 5 orders of 
magnitude which indicates that our gradients are evaluated fairly accurately. Figure 8 
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(a) 



(b) 



FIG. 7. Vorticity fields and streamlines in the numerical solution of (a) direct problem and (b) 
adjoint problem at the first iteration. 

reveals deviations of >Cf {e) from unity for large values of e, which is due to truncation 
errors, and also for very small e, due to round-off errors, both of which are well-known 
effects. These inaccuracies do not affect the optimization process, since the deviations 
observed for very small e are only an artifact of how expression (54) is evaluated, whereas 
large values of e (or, equivalent ly, r) are outside the range of validity of the linearization 
on which the optimization approach is based. We also performed a grid-refinement study 
of the cost functional gradients which indicated that the calculation of the gradients is 
not sensitive to the resolution. Figure 9 shows the decrease of cost functional (47) in the 
case with and without the closure terms a and b in boundary conditions (26e)-(26f) as 
a function of the number of iterations. We observe that the proposed algorithm results 
in a steady convergence despite the complicated nature of the problem, although the 
rate of convergence is relatively slow, especially in the case when the closure model is 
present. 

B. Effective Surfaces for Different Parameters in the Closure Model 

In this Section we employ the computational approach developed in Section V and 
validated in Section VI A to construct effective surfaces corresponding to different values 
of the three parameters {V^, r, T} characterizing the algebraic closure model introduced 
in Section IV. In order to reveal different trends, in Figures 10a, b,c we show the effec- 
tive surfaces obtained by changing one parameter with the other two held fixed. For 
comparison, in these Figures we also include the effective surfaces obtained without any 
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FIG. 8. The diagnostic quantity xv LG (e) as a function of e for different perturbations (★) C(s) = 
sin(27T5/L), (o) C(s) = sin(27rs/L), (*) ('(s) = arcsin(2s/L - 1), (□) C'O) = arccos(2s/L - 1), 
where < s < L. 

io* t , , d 




FIG. 9. Cost functional J(T £i) as a function of the iteration count n for (•) the case without 
closure terms and (■) the case with closure terms a and b in boundary conditions (26e)-(26f). 
The values of the problem parameters are Vd = 1.0, r = 0.25, and T = 3. 

closure model (i.e., with b = in (26f)). The parameters are chosen in such a way that 
the case with {Vd = 1.0, r = 0.25, T = 3} is present in all three Figures 10a, b,c where 
it represents the intermediate solution. First of all, we observe that in all cases smooth 
effective surfaces have been obtained. As regards the results shown in Figures 10a and 
10b, we observe that the effective surfaces approach the effective surface obtained in the 
case with no closure as r — >■ and T — >> oc, respectively. This is consistent with the fact 
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that lim r ^ b = with V d and T fixed, and limy-^ b = with V d and r fixed, cf. (45). 
In the other limits, i.e., for large r (Figure 10a) and small T (Figure 10b), we observe 
that the liquid column (l L becomes much thinner. As regards the dependence on droplet 
velocity V^, from (45) we observe that Vd — ^ < would correspond to the case with a 
vanishing closure model (i.e., b = 0), however, this value of Vd is outside the range of pa- 
rameter values consistent with Assumption 4b. Hence, convergence of effective surfaces 
to the surface corresponding to the case with no closure is not observed in Figure 10c. 
We observe that in the proposed model the closure terms contribute additional flux of 
momentum in the direction tangential to the effective surface which can be interpreted 
as additional shear stress, cf. (26f) and (45). In the cases with large r and small T, 
corresponding to a thinner liquid column Q Ll the effect of the closure model could be 
compared to an increase of the surface tension (although this analogy is rather superfi- 
cial, since the surface tension contributes to the normal stresses). We also add that, in 
addition to the effective surfaces presented in Figures 10a,b,c, for some parameter values 
we also found solutions featuring asymmetric effective surfaces. This nonuniqueness of 
solutions is a consequence of the nonlinearity of the governing system which is reflected 
in the nonconvexity of optimization problem (48) . Since these asymmetric solutions are 
not physically relevant, at least not from the point of view of the actual applications 
of interest to us, we do not discuss them in this work. In problems with multiple solu- 
tions in which such selection cannot be done based on the properties of symmetry, one 
can identify the relevant solution as the one corresponding to the smallest value of cost 
functional J(Tlg) reflecting the smallest residual (46). 

Finally, in Figure lOd, we perform a comparison between the effective free surfaces 
constructed using the algebraic closure model from Section IV and using the time- 
averaged solutions of the original unsteady problem to evaluate the terms a and b via 
relations (23) and (25). The parameters used in this case are Vd — 1.0, r — 0.25 and 
T = 3.0, and the corresponding data for b obtained by averaging over 200 periods of 
droplet impingement in the time-dependent case is shown in Figure 11 (the data for 
the term a is not shown as it vanishes by construction in both cases, cf. (44)). For 
comparison, in Figure 11 we also indicate the values of the components of b obtained 
from expressions (23), and we see that the predictions of the simple algebraic closure 
model developed in Section IV are not too far off from the actual data: as regards 
the vertical component 6 2 , they are within the same order of magnitude (Figure lib), 
whereas for the horizontal component bi the difference is O(10~ 3 ) (Figure 11a). We 
add that for both bi and 6 2 there exists a part of the effective surface T L g, located 
towards the bottom of the liquid column (l L) where the agreement between the closure 
model and the actual data is particularly good. In Figure lOd we note an overall fairly 
good agreement between the effective surfaces obtained with terms a and b evaluated 
in the two different ways. This is, in particular, the case as regards the top part of 
the liquid column Vt L which, somewhat interestingly, does not coincide with the region 
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-' No closure, ' ' r = 0.1, 

•', r = 0.25, ' ', r = 0.5 

with V d = 1.0 and T = 3 
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' ' No closure, ' , V d = 0.5, 

'■■■'V d = 1.0, ' - - - ' V d = 1.25, 
with r = 0.25 and T = 3 

(c) 



— ' No closure, ' ', T = 1, 

T = 3, ' - - - ' T = &, 
with V d = 1.0 and r = 0.25 

(b) 
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' data from the unsteady problem 

(see Figure 11), 
' algebraic model from Section IV 

(d) 



FIG. 10. Dependence of the shape of the effective surface Tlg on the parameters of the 
algebraic closure model introduced in Section IV: (a) droplet radius r, (b) frequency T -1 of 
droplet impingement, and (c) droplet velocity Vd with other parameters held fixed (see captions 
of individual subfigures); Figure (d) shows the comparison of the effective free surfaces Tlg 
obtained using the algebraic closure model from Section IV and the approach described in 
Section VI B in which terms a and b are evaluated based on the actual data obtained from the 
time-dependent flow problem with the parameter values Vd = 1.0, r = 0.25 and T = 3. 
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(a) (b) 

FIG. 11. Comparison of the terms (a) b\ and (b) 62, cf. (25), evaluated based on (solid line) 
the averaged solution of the time-dependent problem and (dashed line) the closure model 
described in Section IV, cf. (45). The data is shown as a function of the normalized arclength 
coordinate s/L along the effective surface Tlg measured from the top, and the values of the 
parameters are Vd = 1.0, r = 0.25 and T = 3. 

where the closure model is the most accurate according to the data from Figure 11. 
We attribute this effect to the nonlinear and nonlocal character of the averaged free- 
boundary problem (26). Finally, we conclude that, despite its simplicity, our closure 
model performs relatively well in the present problem. 

VII. CONCLUSIONS 

In this investigation we revisited the concept of "effective free surfaces" arising 
in the solution of time-averaged hydrodynamic equations in the presence of free 
boundaries 18-20 ' 25 . The novelty of our work consists in formulating the problem such 
that there is a sharp interface separating the two phases in the time-averaged sense, 
an approach which appears preferable from the point of view of a number of possible 
applications. The resulting system of equations is of the free-boundary type and we 
also propose a flexible and efficient numerical method for the solution of this problem 
which is based on the shape-optimization formulation. Subject to some clearly stated 
assumptions, the terms representing the average effect of the boundary fluctuations 
appear in the form of interface boundary conditions, and a simple algebraic model is 
proposed to close these terms (this is to be contrasted with the "classical" Reynolds 
stresses which are defined in the bulk of the flow). 

This work is motivated by applications of optimization and optimal control theory 
to problems involving free surfaces. In such problems dealing with time-dependent 
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governing equations leads to technical difficulties, many of which are mitigated when 
methods of optimization are applied to a steady problem with effective free surfaces. The 
model problem considered in this study concerns impingement of free-falling droplets 
on a liquid with a free surface in two dimensions and is motivated by optimization of 
the mass and momentum transfer phenomena in certain advanced welding processes, 
see Volkov et al. 2 . The computational results shown in this paper are, to the best of 
our knowledge, the first ever presented for a problem of this type where the effective 
boundary has the form of a sharp interface. The computed effective free surfaces exhibit 
a consistent dependence on the problem parameters introduced via the closure model, 
and despite the admitted simplicity of this model, these results match well the effective 
surfaces obtained using data from the solution of the time-dependent problem. 

A key element of the proposed approach is a closure model for the fluctuation terms 
representing the motion of the free surfaces. The model we developed here is a very 
elementary one resulting in simple algebraic relationships. As in the traditional turbu- 
lence research 6 , more advanced and more general closure models can be derived based on 
the PDEs describing the transport of various relevant quantities such as the turbulent 
kinetic energy, the turbulent length scale, etc. In fact, such approaches have already 
been explored in the context of free-boundary problems 20 ' 21 leading to more general, 
albeit less explicit, closure models than the model considered here. 

In addition to investigating such more advanced closure models, our future work will 
focus on quantifying the effect of and, ultimately, weakening the assumptions employed 
to derive the present approach, so that it can be applied to a broader class of problems, 
especially interfacial. At the same time, we will seek to incorporate this approach into the 
optimization-oriented models of complex thermo-fluid phenomena occurring in welding 
processes, such as discussed in Volkov et. al 2 . While the present investigation responded 
to needs arising from a certain class of applications, it has also highlighted a number 
of more fundamental research questions which it will be worthwhile to explore based on 
even simpler model problems such as, e.g., capillary or gravity waves on a flat interface. 
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Appendix A: Characterization of Cost Functional Gradients in Effective 
Surface Calculations 

In this Appendix we obtain expressions for the gradient ^y lg J of the cost functional 
(47) with respect to the position of the effective interface Tlg- Characterization of this 
gradient requires one to differentiate solutions of governing PDEs system (26) with re- 
spect to the shape of the domain on which these solutions are defined. This is done 
properly using tools of the shape-differential calculus 11 ' 12 which are briefly introduced 
below. In the calculations we will assume that the problem parameters {V^,r, T} are 
given. Hereafter primes (') will denote perturbations (shape differentials) of the different 
variables which is consistent with the convention used in the literature 12 . Since fluctu- 
ation variables do not appear in this Appendix, there is no risk of confusion resulting 
from this abuse of notation. 



1. Shape Calculus 

In the shape calculus perturbations of the interface geometry can be represented as 

xfox') =x + t?x / for xef LG (0), (Al) 

where r] is a real parameter, f lg(0) is the original unperturbed boundary and x / : — >• 
IR 2 is a "velocity" field characterizing the perturbation. The Gateaux shape differential 
of a functional such as (46) with respect to the shape of the interface Tlg and computed 
in the direction of perturbation field x' is given by 

JX t „ : x') = lim (A2) 

In the sequel we will need the following fundamental result concerning shape differentiation 
of functionals defined on smooth domains Q(rj, x') and on their boundaries, and involving 
smooth functions / and g as integrands 12 

/ fdn+ [ gda) = [ f'dn+ [ g' da+ 

f + K9 + tt) (x'-n)da, 



Id 



where f and g f are the shape derivatives of / and g ) and k is the curvature of the 
boundary dQ(0). 



2. Differential of Cost Functional 



In order to obtain the gradients ^y lg J of the cost functional (47) with respect to 
the control variable Y L g ) we first need to obtain the Gateaux (directional) derivative 
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of J(Xlg)- Using relation (A3) and substituting / = and g = \\ 2 \ we obtain the 
following expression for the cost functional differential 



J'(T;x') = J r XX'da + J_ Q« x 2 + X (x' • n) da, (A4) 



where 



X ' = n'.K).n + n.(<)'.n + n.(<).n'-n'.(<).n-n.(<)'.n-n.(<}.n' 

r . 9 

- «'7 - 2(p L - p G ) 



(e y -n)(e y • n'). 

(A5) 



Using the following identities of shape calculus 12 n' = -Vr(x' • n) and k! = — A r (x' • n), 
where Vr and A r are the tangential gradient and the Laplace-Beltrami operator, we 
obtain 



J\T-y) = f X ( [n.o-'-n]f-2Vr(x'-n)-K)-n+2Vr(x'-n)-«)-n+ 7 A r (x'-n) 



+ 2(p L - p G ) 



r 



+ { 2 KX +X dn 



(ey -n) (e y • Vr(x'-n)) 



) (x' • n) 



da. (A6) 



Considering Gateaux differential (A6) and invoking the Riesz representation theorem 16 
allows us to extract the cost functional gradient V l 7(f lg) through the following identity 

J , {T;x!) = (vj{f LG ),C) r = / VJCdv, (A7) 
\ / L2(r LG ) Jr LG 

where for simplicity the L 2 inner product was used and ( f = (x! • n) which implies 
that the gradient V J is a scalar-valued function describing the sensitivity to shape 
perturbations in the normal direction. We note that expression (A6) contains terms 
which are already in the Riesz form with the perturbation (x 7 • n) appearing as factor, 
in addition to terms involving the shape derivatives of the state variables, namely (v)' 
and (p) f . The presence of these terms makes it impossible at this stage to use (A6) to 
identify the gradient VJ(Tlg)> In order to transform the remaining part of relation 
(A6) into a form consistent with the Riesz representation (A7) it is necessary to define 
suitable adjoint variables and the corresponding adjoint system. 
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3. Adjoint Equations 

Consider the weak form of system (26) for the variables (v), (v)* G H 1 and (p), (p)* G 



/ [Pl(v) • V(v) - V • «> - p L g] • (v)* - (V • <v»(p)*dx + 

/ \p G (v) ' V(v) - V • «> - PG g] • (v)* - (V • <v»(p)*dx = (A8) 
After integrating the second-order terms by parts, (A8) becomes 

/ [p L (v> • V(v) - p L g].(v)*-(v}-V(p)*+«) : V(v)*dx+ / [p G (v> • V(v) - p G g]-(v) 

-(v).V(pr + K G ):V(vr&-/ n • «> • (v)*^ - / n • «> • (v)*^ 

- / n • (v) (p)*da = (A9) 
Next, using relation (A3) and shape differentiating (A9), we obtain 

/ [p L (v) • V(v)' + p L (v)' • V(v)] • <v>* - V(p>* • (v)' + «>' : V<v)*dx+ (A10) 

/ [p G (v> • V(v)' + p G (v)' • V(v)] • <v>* - V(p>* • (v)' + «>' : V(v)*dx + X = 



where 



i=f j[PL(v}-v(v)- Pi g]-(vr-(v).v(pr+K):v(vr+Kn-«)-(vr 

+ ^ (n • K) • <v)*)+[p G <v> • V(v) - p G gHv)*-<v}-V<p>*+<<T G > : Vv'+k n-«)-(v) 



+^(n.«>.(v>*) ^(x'.n)da- 



n-«)'-(v)*+n'-K)-(v)*+n-«)'.(v)*+n'-«)-(v)^ 



+ n • (v)' (p>* + n' • (v) (p)*J da. (All) 
Performing one more time integration by parts in (A 3) we obtain 

f a | [~Pl(v) ■ V(v>* + p L (v)* • (V(v)) r ]-(v)'-(p)' V.(v)*- /1 (v)'.A(v>'-(v)'-V(p>' J>dx 

+ / ([-p G (v).V(v)* + p G (v)*-(V(v)) T ]-(v)'-(p)'V.(v)*-^(v)'.A(v)* 
Jq g { 

- <v>' • V(p)*Lx + ^ / [(n • V(v>*) • <v>' + (n • (V(v)*f ) • <v>'] da 
J Jr LG 

+ H f [(n • V(v)*) • (v)' + n • ( V(v)*) T • (v)'] da + X = 0, (A12) 
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where (v)* and (p)* can be identified as the adjoint variables with respect to (v) and 
(p), provided they satisfy the following adjoint equations 



p L (v) • V(v)* = p L (v)* • (V(v» T - (V(p>*) - ^A(v)* in fi Lj 



V • (v)* = 



in Q L , 



(A13a) 
(A13b) 



Pg(v) • V(v)* = PG {v)* • (V(v)) T - (V(p>*) - /iA(v)* in ft G , (A14a) 
V • (v)* = in ft G . (A14b) 

Substituting for n' in (All) we can simplify the expression for I as follows 



x(v> • V(v) - p L g] • <v>* - (v) • V(p>* + «> : V(v>* + M . K> • (v) 
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+ — (n • «> • <v>*) + [p G (v> • V(v) - p G g] • <v>* - (v) • V<p>* + «> : V(v) 



+Kn.(cr G ).(v)*+|-(n.(<).(v)*)|(x'.n)^-^ 



[n-(^)'-(v)l2-V r (x'.n).«).(v)* 



V r (x' • n) • «> • <v>* - V r (x' • n) • (v) (p>« 



Imposing the boundary conditions 



<v>* r = (v)* 



<v>* = 0, 

the expression for the differential of the cost functional becomes 



da. (A15) 



(A16a) 
(Al6b) 



J'(T, x') = / i - V r (x' • n) • «) • n - V r (x' • n) • «) • n + 7 A r (x' • n) 



+ 2(p L - p G ) 



(e y • n) (e y • V r (x' • n)) 



[ Pi (v) • V(v) - p L g] ■ <v)* - (v) • V(p>* + «) : Vv* 



d 



+,n.K).(v)*+-(n.«)-(v)*)+[p G (v) • V(v) - p g]-<v)*-<v>.V<p>*+<<r£> : V<v)< 



+ K n> G )>)* + |-(n. <<)>>*) 



(x'-n)U<7 (A17) 



which is now consistent with Riesz's representation (A7). Finally, after applying tan- 
gential Green's formula 12 to the terms involving Vr(x' • n), the cost functional gradient 
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can be identified as follows 
V JiT LG ) = 



+2(p L -p G ) 



- n- V r • «) - n- V r • <<r G ) + Qk X 2 + X g ) 

(e y -n)-[p L (v) • V(v) - p L g]-(v) 



^2 



(v).V(p)*+«) : V(v}*+ K n-«)-(v)*+^ (n • «) • <v)*)+[p G (v) • V(v) - p G g]-(v)^ 



(v) • V(p>* + «) : V(v)* + k n • «) • <v>* + ^ (n • <<r G ) • <v>*) +7^ 



on r LG - 

(A18) 



Consistency of this expression for the cost functional gradient is demonstrated in Section 
VIA. 
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